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ABSTRACT 



Large matrix storage constitutes a limitation on the 
applicability of most numerical techniques including the 
Finite Element Method, when very accurate results are required. 
This is particularly true when dealing with Boundary Value 
Problems. In order to surpass this difficulty a new method 
to solve these problems has been devised which does not re- 
quire matrix storage while still providing the possibility of 
accuracy improvement. 

Although restricted to one-dimensional, linear differen- 
tial equations of the form Y v ; (x) = f(x) this new approximating 
technique gives acceptable results. The method will perform 
equally well for problems with exact or non-exact integrable 
forcing functions, continuous or discontinuous, or functions 
existing only as a set of values at discrete points. 
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INTRODUCTION 



I . 

For many engineering problems it is not always possible 
to find an exact solution. An exact solution is an analyti- 
cal mathematical expression that gives the value of the unknown 
at any point within a previously specified range. For problems 
involving complex material properties and boundary conditions, 
different numerical techniques have been developed to approxi- 
mate the exact solution to a more or less acceptable degree 
of accuracy. One of these techniques is the Finite Element 
Method, where a given problem is to be discretized, that is to 
say, approximate solutions are to be found at discrete points 
in the body. The way one accomplishes this, is by subdividing 
the whole body into finite elements. The solution is then 
formulated for each small unit and combined to obtain the 
solution of the whole system. Obviously the greater the number 
of elements, the better the accuracy. 

However, despite the fact that high-speed digital computers 
have enabled engineers to successfully apply these numerical 
techniques, we still face the problem of large matrix storage. 
Complex problems requiring very accurate approximate solutions 
lead to large matrices and consequently larger computer memory 
will be required. This fact constitutes a limitation on the 
applicability of these numerical techniques. 
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The method we are dealing with in this work intends to 
solve this difficulty by reducing a B.V.P. to an I.V.P. Unlike 
the "shooting" method, where basically the same idea is used 
together with an iterative scheme to achieve a solution, the 
method developed here achieves a solution without iteration. 

Although restricted to a particular type of linear differ- 
ential equations and only to one-dimensional problems, as it 
will be shown, this method can be applied to many problems 
where no exact integrable forcing functions occur, or, even 
more, the function exists only as a set of values at specific 
points . 

Basically what the method does is approximate the exact 
solution by a set of linear functions, each of them applying 
in a very small interval. However, these functions are not 
independent from each other. We use the previous one to find 
the next, until we eventually reach the other end of the range 
in consideration. Depending on the order of the differential 
equation, we will find as many sets of linear functions as 
required by its order, each set applying to a specific inte- 
gration, and as before, every set depends on the previous one. 
The way in which we will do this will allow us to correlate 
all possible boundary conditions in an explicit manner which 
transforms boundary value problems into initial value problems 
and by doing that, get to the result. 

The method is indeed an approximation and as such, it is 
subject to error. We can say, however, that by decreasing the 
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step size we get better accuracy, provided we carry enough 
significant digits to take care of round-off errors. With 
respect to the later a good programming technique is required. 
Finally we shall apply the method to several practical applica- 
tions, specifically to beam problems where a fourth order 
differential equation occurs, and a variety of boundary 
conditions can be given. 
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II. BACKGROUND 



A. GENERAL 

In this section we are going to introduce the basic ideas 
of the method we will be dealing with. Let 

Y’(x) = f(x) (2.1) 



where it is understood that the function f (x) may or may not 
be exactly integrable, but we can integrate it numerically. 
From the above relationship we have: 

Y(x) = / f (x)dx + Y q (2.2) 

Here the constant of integration Y q is assumed to be zero 
for the moment. In the most general case, it is clear that 
there is no way to know what the function Y(x) is, since we 
may be dealing with non-integrable functions. However, we 
can approximate the function Y(x) by a linear function y(x), 
provided the interval in which this approximation applies is 
sufficiently small. So we take 

Y ( x) = y(x) = mx + c (2.3) 

Now, by following the general procedure of integration, and 
from (2.2), we have: 
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(2.4) 



b b 

mx + c | = / f ( x ) dx 

a a 



where (a,b) denote the limits of the interval under considera- 
tion. After simplification it is found that: 



b 

/ f ( x ) dx 



(2.5) 



that is to say, the slope of the approximating line given by 
(2.3) can be explicitly determined. Actually, as we will see 
later, it can be said that the "general" approximate solution 
to (2.1) has been found within the interval from a to b . In 
order to determine the "particular" solution, in other words 
to find the intercept c of (2.3), an auxiliary condition must 
be imposed. Let 



Y(a) = Y q 



then : 



Y(a) = Y q = ma + c 



or 



y(x) = m(x-a) + Y Q a <_ x <_ b 



( 2 . 6 ) 
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where m is given by (2.4) . Now we have fully determined the 
approximate solution inside the given interval. Note that 
(2.6) will give results at points a and b, as accurate as the 
numerical integration performed on (2.5). 

Fig. 2.1 shows the whole process so far. The upper curve 
represents f(x), which is given. The lower curve is Y(x), 
the integral of f(x); however within the interval a to b, Y(x) 
is approximated by a line. Point Y q is the given auxiliary 
condition and point Y^ can be determined exactly from (2.6) , 
by letting x = b. 

The next step now becomes evident, since point Y^ can be 
determined. We are now in the same position as before, so 
all we need to do is repeat the process over again. However, 
this time with a new auxiliary condition; the last one we 
have just found, and over the next interval. We keep going 
this way until we reach the other extreme of the range where 
the differential equation applies. See Fig. 2.2. 

In summary, we can say that every step we take we are 
solving a "new" differential equation by approximating lines 
which gives us true values at the points of intersection. 

The fact that we have a "new" differential equation at every 
single step allows us to deal with discontinuous functions and/ 
or functions existing only as a set of values. 

B. THE FIRST INTEGRATION 

At this point let us introduce a new parameter: 

h = b-a (2.7) 
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FIGURE 2.2 
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where h is assumed to be small and constant throughout the 
analysis. We are now in a position to determine a general 
relationship which will allow us to find the value of Y(x) at 
discrete points , namely at the intersection of the straight 
lines. From (2.6) we have: 

y(b) = m(b-a) + Y = mh + Y 
1 o o 



In general. 



Y. . = m.h + Y. . 0 < i < n 

l-l l l-l — — 



( 2 . 8 ) 



where 



I f ( x ) dx 



m . 



k i-l 



(2.9) 



but : 



Y. . = m. .h + Y. _ 

l-l l-l 1-2 



Y i-2 = m i-2 h + Y i-3 



Y. = m.h + Y 
1 1 o 



Replacing these last relationships in (2.8) and after factor- 
ing h, we have: 
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or 



Y r> = (m, + m~ + ... + m )h + Y 

n 1. z, no 



n 



n 

h Y m. + Y 
i=l 1 



( 2 . 10 ) 



Similarly, 



m-, 



m. 



h 



/ f(x)dx 



2h 

/ f (x) dx 



Adding terms gives:.’ 



n 

l 

i=l 



m . 



0 



nh 

/ f (x) dx 



( 2 . 11 ) 



We can replace in (2.10) to obtain: 



n 

Y = / f ( x ) dx + Y 

n 0 ° 



( 2 . 12 ) 



since x = nh . Equations (2.10) and (2.12) constitute two very 
important results. Eq . (2.12) is indeed not an unexpected one, 
and several considerations can be drawn from it. First of all 
the value of Y n does not depend on the step size h, at least 
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for this first integration, so its accuracy is determined only 
by the exactness of the integration routine used. Second, the 
values of Y(x) can be directly determine at any arbitrary point 
n, that is to say we need not find its previous values as is 
the case for many numerical methods. Later this fact will 
become more useful. Finally, we see that, by transposing terms 
we can find either the left or right B.C., and this apparently 
unimportant fact is indeed a key step, since as we will see 
later, we will be able to correlate B.C.'s of higher order 
differential equations and solve for them without having to 
generate all previous values of the unknown function. In other 
words, we could transform a Boundary Value Problem into an 
Initial Value Problem and vice versa, depending on what is 
known, and what we are looking for. 
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Ill . HIGHER ORDER DIFFERENTIAL EQUATIONS 



A. GENERAL 

The previous chapter was dealing basically with the main 
ideas of the method, and we have solved a first order differ- 
ential equation. We are now going to extend the method to 
higher order equations. Essentially what we will do is, inte- 
grate several times the right hand side of the given equation; 
as required for its degree, using the same approach as before. 
However, we should keep in mind that it is only for the first 
integration that we will use the given function f(x). After 
this integration, this function is no longer available because 
what we have is a set of straight lines, each applying to a 
specific interval. We need to find now a method that will 
allow us to perform a second integration of the function f(x), 
from the given set of lines. 

In order to do this, let us consider the curves of Fig. 
3.1. Suppose we have been given the following second order 
differential equation: 



Y" (x) = f 2 ( x ) 



(3.1) 



represented in Fig. 3.1 by the upper curve. Let us assume we 
have determined Y'(x) = f-^(x) after a first integration by the 
method just introduced, Let 
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FIGURE 3.1 
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Y(x) 



f (x) 



be the exact result, which is the lower curve in the same 
figure. Now let 



y ' (x) = m-^x + c-^ 



and 



y ( x ) = mx + c 

be the approximating lines to f^(x) and f(x), respectively, 
that applies inside the interval from a to b only. In this 
figure, we can easily see that the area enclosed by the upper 
straight line is almost equal to the area enclosed by the 
function f^(x) . Obviously as the interval becomes smaller and 
smaller, both areas tend to be equal; in the limit they are 
indeed equal. It is evident now that we can approximate the 
exact integration of f-^(x), by integrating the line 
y' (x) = n^x + c. , provided the interval is small enough so the 
error is negligible. Now, as before, we need to determine m 
and c. Recall that m^ and c^ have already been found. So by 
virtue of (2.5): 

b 

/ f ^ (x) dx 
a 



19 



but by (2.4) 



b 

/ (m^x + c^)dx 



or : 



1 

m = —(a + b) + c 1 (3.2) 

So we know now the slope of the next line; the line belong- 
ing to the solution function. In order to find the intercept 
c, another auxiliary condition must be supplied. Let: 

Y(a) = Y a 

then, replacing this condition in the line y(x) = mx + c , we 
get : 



Y(a) = Y = ma + c 

d 

It can be shown that: 

y(x) = m(x-a)+Y (3.3) 

3 . 

where m is given by (3.2). Eq . (3.3) gives the approximate 

solution of the original second order differential equation, 
inside the interval a to b. Note that both auxiliary conditions 
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have been given at the left hand side. However, it really 
need not be like this, as will be shown later. 

At this point some considerations are in order. With 
reference to Fig. 3.1 where the upper curve is the given 
function Y"(x) = f 2 (x) , the middle curve is Y'(x) = f 1 (x), 
and the lower one is the solution Y(x) = f(x) . 

As explained previously we have approximated the solution 
by straight lines. Consider now the slope of the lower 
straight line m which is given by (3.2) . This equation can be 
written as : 

m = ■^•[(m^a + c^) + (m^b + c^) ] 

m = i ty 1 (a) + y^ (b) ] 

m = l tY a + V (3.4) 

Since the ordinate at any point x in the Y'(x) curve is 
indeed the slope of the Y(x) curve at the same point, then 
by virtue of (3.4) , we can say that the slope m of the lower 
straight line is the average of the slopes corresponding to 
its end points, namely point Y q and Y^ in Fig. 3.1. Recall 
that these two slopes are given exactly by the ordinates of 
the middle curve Y'(x). Furthermore let us evaluate 
y'(x) = m-^x + c-^, the equation of the upper straight line, at 
the midpoint of the interval, that is we let: 
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X 



^-(a + b) 



in the line y' (x) and we get: 



y' ( 



a + b 



) = 



( 



a + b. 



+ c. 



= m 



As seen, it reproduces m, the slope of the lower straight 
line . 

Now we can have a deeper insight of the whole process, 
and why it works. See that what we are really going to do is 
to use the slope at the midpoint (the average of the slopes at 
the extremes) to approximate the true slope of the straight 
line y(x), and later determine the corresponding intercepts. 

In fact, had we known exactly the function Y'(x) = f-^(x) , we 
would have been able to determine the true slope of y(x) by 
direct integration as before. As we can see now, we are intro- 
ducing a source of error, and from now on we will not have as 
accurate results as in the first integration. 

B. DEVELOPMENT OF THE METHOD 

Before we go into a general development of this method a 
slightly different notation must be introduced, and we will 
be using it throughout the rest of this research. The type 
of differential equations we will be looking at are of the 
form: 



Y 1 V (x) = f ( x ) 



( 3 . 5 ) 
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Such equations arise in the study of deflections of beams, 
and we shall apply the method to practical beam problems. Let 



us introduce the notation to be followed. We will write the 
approximate solution line of ( 3 . 5 ) as: 

y ( x ) = mx + c 

in a clean way without subscripts. The approximation line 
to the first derivative will then be: 

y ’ (x) = m^x + c^ 

and for the second derivative: 

y"(x) = ir^x + C2 



and so on. 

The boundary conditions will be represented as: 



Y (Xq) 


= Y 

0 


Y(x n> 


= Y 

n 


Y'Uo) 


= Y’ 
0 


Y’K) 


= Y’ 
n 


Y " (x ) 
0 


= Y" 
0 


Y" (x n ) 


= Y" 
n 



etc., where the subscript o refers to the left hand side, and 
the subscript n to the right, of the range under consideration. 
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There is another subscript that we must keep track of. 

Recall that the whole solution function is going to be repre- 
sented by a set of approximating lines, see Fig. 2.2. Then if 
we start at the left end, we shall use the subscript n as a 
second indice to identify each line as we move rightwards. Be 
aware that we will take the subscript i immediately to the 
right of each line to name it locally. 

Some words must be said about the auxiliary conditions. 
Recall that in a specific problem the number of these condi- 
tions are equal to the order of the differential equation to 
be solved. If all these conditions are given at the same 
point, then we are dealing with an Initial Value Problem. 
However, if the auxiliary conditions are given at both ends 
of the range of interest, then we have a Boundary Value Prob- 
lem. When we are able to have an exact solution, we usually 
are led to a set of simultaneous equations involving the unknown 
B.C.'s; those can be solved algebraically and lead to the 
whole solution. However for equations with no exact solution, 
like the ones we are interested in, there is no way to deter- 
mine the conditions at the other end of the interval, and that 
is precisely the goal of our approach. We need to correlate 
in an explicit way the known conditions and the unknown ones 
by means of equations of the form of (2.10) , so we can be able 
to find them. 

n 

Note that in this equation the term h £ m. is independent 

i=l 1 

of B.C.'s and it can be readily found even if we do not know 
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them. By similarity to an exact problem we can look at this 
term as being the general solution of the differential equation. 
Note that this term actually links the left and the right 
boundary conditions in a very explicit way. This is the 
heart of what the present method is all about. With this 
idea in mind we are now going to find several relationships 
for successive integrations. 

Let us start by noticing that (2.10) is a general relation 
that applies not only to the first integration but also to 
successive ones. The only term that is particular for each 
integration is the summation of the slopes. So in order to 
perform the next integration we shall determine this term 
specifically and substitute it into (2.10). We are now dealing 
with a fourth order differential equation since this type 
occurs in beam problems. The application to other orders will 
be self-evident. Let 



Y V (x) = f 4 (x) (3.6) 

After a first integration we have: 

n 

Y"' = h y m. + Y"' (3.7) 

n l 3 i o 

i=l 

where the subscript 3 in the summation indicates that these 
slopes belong to Y"' = f^ (x) and are given by (2.11). We need 

to find now 
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(3.8) 



Y 



II 

n 



n 



h 1 2 

i=l 



m . 



1 



+ Y" 
o 



Following our convention, m 9 . is given by: 

z , ± 



m 



m 



2,i 



. , + x . ) + b, . 

2 l-l l 3,i 



but 



)^ . = Y'." , - m 0 . x . , 

3,i l-l 3,1 l-l 



then 



m 



m 



2,i 



3 , i 



(x, 



i-1 



+ x . ) + Y V 



i-1 



- m 



. x . 



3,i i-1 



or 



m 



2 , i 




+ 



v " 1 
1-1 



(3.9) 



Now, since 



v M 1 
1 n-l 



n-1 

h 

i=l 



m . 

l 



+ Y"' 
o 



then 
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(3.10) 



m 



2 , n 




n-1 

+ h y, m. + Y'" 
L 3 1 o 

i=l 



This last relationship in the way it has been derived 
becomes a general one, and expresses the slope m 2 n at any 
point i as a function of the previous slopes. Similar rela- 
tions can be derived for successive integrations just by 
shifting the first indice. Now we need to find the summation 
of these slopes. It becomes: 

n , n n-1 i 

y„ m. = % V, m. + h y Y m. + nY"' (3.11) 

L 2 1 2 L 3 1 L 3 3 o 

i=l i=l 1 j=l 

Now, replacing (3.11) in (3.8): 



Y" 

n 



l n-1 i 

t x 4 y Q m. + h y y m. + nY ' M ] + Y " 

2 ^3 1 -l-i L 3 3 o 0 

i=l j=l 



htt 



or 



Y" 

n 



1 n n-1 

h [ l ^3 m i + .1. 

i=i 1=1 



y o m.] + nhY + Y" (3.12) 

L 3 1 o o 



j = l 



Let : 



1 n 
= k l 



n n-1 i 

3 m i + l 



m . 

3 3 



i=l i=l j=l 
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Note that, as stated above, this last term S 2 is also 
independent of boundary conditions and can be readily deter- 
mined. Equation (3.12) gives the particular solution to a 
second order differential equation as well as, in this 
specific case, it represents the second derivative of Y(x) . 

The whole process, so far, can be summarized as follows: 
Equation (2.10) gives the solution at any point n of any 
integration, first second, etc. In order to use this equation 
we need to find the summation term, and this is given by a 
relation of the form of (3.11) , which is also a general rela- 
tionship and expresses the summation of slopes in terms of 
the previous one . Recall that the only slope summation we know 
is the one involving the known function f(x) , and is given by 
(2.11). That is why all further summations must be expressed 
as a function of this one. This last summation is then re- 
placed in (3.8), and simplified if possible. The next step is 
to identify the terms that are independent of boundary condi- 
tions and isolate them as in (3.12) and to evalute them 
separately . 

We can continue in the same way until we get a solution 
formula for Y(x). It can be shown that the complete solution 
for a fourth order differential equation of the form Y 1V (x) = f(x) 
is given by the following relationships: 





(3.13) 




hS„ + nhY 

A 



+ Y 



O 



II 



(3.14) 



o 



28 



(3.15) 



Y' = h 2 S, + N 1 h 2 Y"' + nhY " + Y' 

n 1 1 o o o 



'n 



h 3 S + N h 3 Y"' + N,h 2 Y" + nhY' + Y 

o o o 1 o o o 



(3.16) 



where : 



n 

S, = J m. 
3 i=i 1 



(3.17) 



n 



■i n-1 i 

= \ l m i + l l m -j 

i=l 1 i=l j=l 3 



(3.18) 



•. n n-1 j n-2 i j 

i I m i + I t **+ l l I 

i=l i=l j=l 3 i=l j=l k=l 



m, 



(3.19) 



■. n 0 n-1 i 0 n-2 l i n-j i j 

; = i l + t l l + t l l 1 + I J I 

i=l i=l ]=1 J i=l j=l k=l i=l 3=1 k=l 1 _ 



n-2 i 



(3.20) 



n-1 

N x - I i 

1 i=l 



(3.21) 



n-1 n-2 i 

N 0 = £ + l i + I. I j 



i=l i=l j 



(3.22) 



and 
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II 



1 

h 



ih 

/ f(x)dx 

(i-l)h 



0 < i < n 
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IV. ERROR ANALYSIS AND CONVERGENCE 



Refer to Fig. 4.1 where we have plotted the functions 
Y' = f ^ (x) and Y = f (x) together with their respective 
approximating straight lines inside the interval x = a to 
x = b. We assume that the upper straight line belongs to a 
first integration; in other words we have not introduced any 
error so far. A second integration, however, will be performed 
using this line y-^(x) since we do not know the function f^(x). 
It is worth considering now what happens when we integrate 
the line instead of the function itself. Let m be the slope 
of the lower straight line. This is given by; 

b 

/ f 1 (x)dx 



Let m* be the approximating slope given by: 



/ (m^x + 0 -^) 



dx 



m 



* — 



b-a 



(4.1) 



Now for the particular case shown in the figure, we have; 



b b 

/ (m^x + c^)dx > / f^(x)dx (4.2) 

cl 3. 
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FIGURE 4.1 



32 



since the area enclosed by the curve is smaller than the area 
enclosed by the line (we assume that there are no inflection 
points and/or discontinuities inside the interval). It is 
now clear that: 



m* > m 



and the error introduced is: 

b b 

e = / (m.x +c.)dx - / f. (x)dx (4.3) 

a-b a 1 1 a i 

which is equal to the shaded area in Fig. 4.1. 

As indicated in Chapter III, m* comes out to be the average 
of the slopes at the end points of the lower straight line, 
and it is this slope we work with. So, the actual approximat- 
ing line y*(x) = m*x + c*, shown in Fig. 4.2 is steeper than 
y(x) since m* > m. The distance ED becomes the error introduced 
and we can not evaluate it because we do not know f(x). Another 
integration using y* (x) will make things even worse. Note 
that the error to be committed this time will be larger and 
is given by the shaded area in Fig. 4.2. 

So far it appears that the method is divergent in nature 
due to the fact that the error will grow bigger and bigger 
unless we find some means to correct it. At this point, the 
only possible way to do this is by reducing the "shaded" 
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FIGURE 4.2 
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areas, which means that we have to reduce the step size. 
Fortunately, we can say, at least in theory, that in the limit 
there will be no error; in practice however, a very small 
step size will increase the round-off errors. The examples 
provided show that the approximation is indeed acceptable. 

There is a special case, however, where an error correc- 
tion can be made. This applies only to boundary value problems. 
Recall that a solution of any differential equation by this 
method is given by the general relationship: 



Y 

n 



h 



n 

l 

i=l 



m . + 
1 



In B.V.P., we know the conditions at both ends of the interval. 
So this last equation can be written as: 



n 



l 

i=l 



m . 



1 




( 4 . 4 ) 



where the right hand side is known. Now since we are actually 
dealing with approximations to the true slopes, what we really 
have as a summation is the next term, call it M: 



n 

M = 7 m* 

i=l 1 
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Our approximate solution is then: 



n 

Y* = h J m* + Y (4.5) 

n . L , o ' ' 

1=1 

and the error at point n will be given by: 

e n ’ Y S - Y n (4 ' 6) 

Now let us see how the correction could be carried out 
individually at each interval. To do this we must find the 
exact slope m. It is given by: 

m = m* - e 

where e is the shaded area of Fig. 4.1 and is given by (4.3). 
The slope summation is therefore: 





je 



Multiplying by h and adding Y q to both sides, we obtain: 



h 



i=l 



m. + 

l 



h 




+ Y 

o 



jhe 



or 
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Y . 

3 



Yj - jhe 



(4.7) 



When j = n, we have: 



Y = Y* - nhe 
n n 



or 



Y* - Y = nhe 
n n 



but from (4.6) we get: 



e 




Replacing this last relation in (4.7): 



(4.8) 



Y . 

3 




i 

n 



e 



n 



and the general relationship given by (2.10) transforms into 



Y . 

3 



3 



h 7 m? 
i=l 1 



+ Y 

o 




(4.9) 



which is the corrected solution at any point j. In summary, 
what we need to do to correct the solution is: a) determine 
Y* by performing the integration without any correction by 
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using equation (4.5); b) compare Y* with the given B.C. Y n 
and obtain e; c) Perform again the integration this time using 
the relationship given by (4.9). 

From (4.9) we can find a whole set of corrected equations 
similar to those given by (3.13) to (3.20). However, by doing 
this we will enormously complicate those relationships and 
the computational effort, together with round-off errors, may 
not give any advantage at all, especially since we expect the 
local error e to be very small. 

There is, however, another stronger reason not to do that. 
The fact is that we have assumed that the error e is a con- 
stant and this is not always the case. If we restrict our- 
selves to integrate linear and/or constant functions, then the 
corrected method could be justified since for these functions 
the error given by (4.3) is a constant, for the shaded area 
of Fig. 4.1 will always be the same. In more general cases, 
the error e will be unpredictable and no correction can be 
made . 

This example serves only to illustrate how the error 
correction could be carried out. But it would not be practical 
and a reduction in step size appears to be the best correction 
as we will see in the examples. 

As stated above, the method appears to be divergent in 
nature, so our approximate solution will look like Fig. 4.3. 

In this figure, note that for any integration after the first 
one, the set of approximating lines diverge from the exact 
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solution, the further we move to the right, the greater the 

error introduced. In dealing with I. V. P.'s, the points Y' 

o 

and Y q are known and these are our starting conditions. The 
solution will appear then as shown in this figure. Since we 
do not know the exact end points at the right hand side, 
namely points Y' and Y, there remains an uncertainty in the 
accuracy of the solution. 

For B.V.P.'s, however, the other extreme points are known. 
Referring to Fig. 4.3, suppose we are given as boundary condi- 
tions, points Y^ and Y q in the lower curve of this figure. 

The method we are dealing with requires that we know the 
initial points to be able to start the integrating algorithm. 

If we know the exact starting points, namely Y^ and Y q , then 
by using these initial conditions, our solution will appear 
as shown in Fig. 4.3. That is to say we will not end up at 
point Y which is exact, but at point Y*. To be able to reach 
the exact point Y, which we assume is known, we need to give 

the algorithm a "wrong" starting point Y'*. It is this approxi- 

o 

mate starting point which we find using the relationships given 
by equations (3.13) to (3.22) that allow us to match the exact 
solution in the whole range as shown in Fig. 4.4 in the lower 
curve . 

Note that if we were using corrected relationships we 
would be able to supply the algorithm with "exact" starting 
values, and still match the exact end points. Unfortunately, 
it is now the upper set of solution lines that absorb the 
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inherent error of the method since we have given them approxi- 
mate initial values. The solution now is going to look like 
Fig. 4.4. The lower solution becomes "exact" but the upper 
actually shifts apart from the exact. We can explain this by 
saying that the known B.C. 's are fixed by the algorithm while 
the free (unknown) ones will take the error, and this is 
exactly what happens in the examples in the next section. In 
summary, all we do is to shift the error from one place to 
another. We can not get rid of it unless we use corrected 
relationships or an infinite number of intervals to minimize 
the local error. 
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V. APPLICATIONS AND RESULTS 



A. GENERAL 

In this chapter we are going to study a particular appli- 
cation of the integrating algorithm introduced in the previous 
sections, namely, the deflection of beams, where a fourth 
order linear differential equation occurs, and several combina- 
tions of B.C.'s can be considered. 

From Mechanics of Solids, the relationships governing the 
deflection of beams are given by: 

i v 

EIY (x) = q(x) (load intensity) (5.1) 

EIY"' (x) = V(x) (shear force) (5.2) 

EIY" (x) = M (x) (moment) (5.3) 

Y ' (x) = slope 

y(x) = deflection 

The method we are dealing with is basically a method of 
successive integration. Because of that we can start to inte- 
grate from any of the above relationships provided we know 
explicitly the right hand side of the equation. In some cases. 
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a given problem could be solved by integrating the moment 
second order differential equation. By doing this, we would 
be able to achieve more accuracy since only two integrations 
need to be performed and the accumulated error will be rela- 
tively small . 



However , 


when the loading of the beam is a complicated 


dis tribution 


and the expression for the bending moment is 


difficult to 


obtain, then the fourth order differential equa- 


tion will be 


the one to use. One advantage in using this 



expression is that the integrating algorithm will provide 
results for the shear force, bending moment, slope and 
deflection simultaneously. When using the second order equa- 



tion we only 


get slope and deflection. 


Since eq, 


. (5.1) is a more general relationship, we are 


going to use 


this expression in the examples. In fact, most 


problems can 


be expressed in that way. Recall that a concen- 


trated force 


can be treated as a distributed force acting in 


a very small 


interval, and a concentrated moment reduces to a 



couple of concentrated loads acting in opposite directions a 
small distance apart. Example 2 illustrates this procedure. 
The examples provided in this chapter account for all 



types of B.C, 


. 's. Fig. 5.1 shows the 4 cases to be treated. 


and we shall 


study them one by one in the next subsections. 


In all cases 


the flexural rigidity El and the length L is 


set equal to 


1, for simplicity. Discontinuous types of loading 



are emphasized, and of course, many combinations of loadings 
can be dealt with by simply using the principle of superposition 
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FIGURE 5.' ^ 
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The computer program used to solve these problems is shown 
at the end of this chapter, and it is a straightforward coding 
of equations (3.13) to (3.22) . This programming technique, in 
particular, tries to minimize round-off errors. For now, it 
suffices to say that all problems are solved using SINGLE 
PRECISION, one hundred intervals (problem 2 uses 200), and we 
get accuracy to the fourth and even to the fifth decimal places 
in some examples. Finally, only the fourth integration, the 
deflection, is compared with the exact solution. Recall that 
the error is cumulative, and it is here that we expect the 
bigger error. We may conclude then that preceding integrations 
are more exact. However, we should keep in mind that an 
approximated starting point is needed in some cases, and the 
results will be shifted by some amount from the exact. We 
will see this as we proceed into the next section. 

B . EXAMPLES 

1 . Example 1: Clamped-Free Beam 

Refer to Fig. 5.1a. This is a statically determinate 
cantilever beam loaded over half its length by a uniformly 
distributed load of 1 unit of weight per unit of length as 
shown. We shall determine the shearing force, bending moment, 
slope and deflection at discrete points of the beam. 

2 . Solution 

This is a B.V.P. but for purpose of illustrating how 
to use the method when dealing with I. V. P.'s, we will treat 
it as such. All the required initial conditions can be drawn 
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from the given configuration. Here the function f(x) to be 
integrated is discontinuous and given by a set of two differ- 
ent functions defined by: 

f 2 ( x ) =0 0 x < 0.5 

f 2 ( x ) = 1.0 0.5 _< x £ 1.0 

We need to determine the initial values to start the 
algorithm. From statics, it can easily be shown that the 
initial conditions are as follows: 



Shear : 


Y»» — 

o 


-0.5 


Y 1 = 

n 


o 

• 

o 


Moment : 


Y" = 
o 


0.375 


Y" = 
n 


o 

• 

o 


Slope : 


Y’ = 
o 


0.0 


Y 1 = 

n 


Unknown 


Deflection : 


II 

0 


o 

• 

o 


Y 

n 


Unknown 



By supplying these initial values to the program, we 
obtain the computer output shown on the next page. The last 
column shows the exact deflections. In these results, note 
that all the initial values are exact, and that the last value 
found, namely the deflection of the beam at x = 1.0, is expected 
to have the biggest cumulative error. However, as we can 
see, we get 3 to 4 digits accuracy. Here we can conclude 
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EXAMPLE 1 . 

CLAMPEC-FREc SEAM 



BOUNDARY 


CON C I T IC.N S : 


Y ( 0 ) = 


O.CCOO 


Y (C)= 0.0000 








M(0 ) = 


0.3750 


V(0 )= -0.50C0 




x : 


SHEAR 


MOMENT 


SLOPE 


C SELECTION 


EXACT OEF 


0. JO 


-0.5CCC00 


0.375000 


O.COOOOO 


C. 000000 


-0.000000 


0.05 


-0.500000 


C. 350000 


0 .018125 


C. 000453 


C. 000458 


0.10 


-0. 500000 


C. 325000 


U.C35000 


C. 001792 


0.001792 


0.15 


-0.500000 


C. 300000 


0 .050625 


C. 003937 


0.003938 


0.2 G 


-0 . 5GC0Q0 


0.275000 


0 .065000 


C. 006333 


0.006333 


0.25 


— C . 5 C C COO 


0.2 5 OCO C 


C.C73125 


C. 010416 


0.010417 


0 . 3 C 


-C. 5 C C COO 


0 .225000 


0.090000 


C. 014625 


0.014625 


0.35 


-0.500000 


0.200000 


0. 100625 


C. 019395 


0.019396 


0.40 


-0. 5CCC00 


0.175000 


O.llOuOO 


C. 02466o 


0.024667 


0.45 


-0.500000 


C . 1 5u00 0 


0.118125 


C. 030375 


0. 030375 


0.50 


-0.50 0000 


C. 125000 


0. 123000 


C. 036453 


0.036459 


0.5 5 


-0.45 CC JO 


0.101250 


0 • 130646 


C .04285-+ 


0.042355 


0 . b 0 


— 0.40 COOO 


C .C8U000 


0.135167 


C. 0+9504 


0.049505 


0.65 


-0.3 5CC0C 


0.061250 


0 .138638 


C. 05 6354 


0 .35o 355 


0 . 7 C 


- 0.3 0 COO 0 


0 . 04>00 0 


0 .141334 


C. 063353 


0.063 358 


0. 75 


-O.25CGO0 


0.031250 


0. 143230 


C. 070475 


0.070475 


0 .6 C 


-0. 2CCC00 


C .020000 


0 . 144501 


C. J77o70 


0 . 077 o7 1 


0.25 


-0 .150000 


0.011250 


0.145272 


C. 034916 


0.034916 


0.90 


-0. lOCCuO 


C . C 05000 


0 . 1456o7 


C. 092191 


0.092 191 


0.95 


-0. C5CCOO 


u .00 1249 


0. 145313 


C. 099479 


0.099479 


1 .00 


-0 .000000 


-C. 00000 1 


0 . 145334 


C. 106770 


0.106 770 
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that previous values of all other columns are even more 
accurate . 

3 . Example 2; Clamped- Roller Beam 

Refer to Fig. 5.1b. This is a statically indeterminate 
beam loaded by a moment of 10 units of weight x unit of length, 
acting at the middle of the beam. As before we are going to 
determine the shear, moment, slope and deflection of the beam. 

4 . Solution 

This example will emphasize three things: How to deal 

with a) concentrated forces, b) concentrated moments and c) the 
necessary initial conditions required for this particular 
case . 

a. Concentrated Forces 

We can conceive of a concentrated force as a dis- 
tributed force of high intensity distributed over a very small 
length of beam. Recall from the previous chapters that the 
slope of any approximating line at any point i is given by 
(2.9) where f(x) is the distributed load for the case of beams. 
In this equation the integral in the numerator represents the 
area under the f(x) curve. This area is actually the total 
load in the specified interval. By decreasing the interval 
while still holding the same inside area constant (same total 
load) , the distributed force tends to a concentrated force of 
equal magnitude as the total load. In the limit it becomes a 
concentrated force acting on a mathematical point on the beam. 
For our purpose, in the integrating algorithm we simply let: 
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0 



for 



m . 

1 



x ¥■ x . 
1 



m . 

l 



P 

h 



for 



x = x . 

l 



/ f(x) dx = P for x = x. 



where P is the concentrated force acting at point x^. 
b. Concentrated Moments 

A concentrated moment can be treated in a similar 
way if we use the second order differential equation. However, 
in our case we need to decompose the given moment into a pair 
of concentrated forces of equal magnitude acting in opposite 
directions at equidistant points from the location of the 
given moment. From statics we know that: 



Moment = Force x Distance 



Here we can have several combinations of force 
times distance provided we keep this product a constant equal 
to the given concentrated moment. However as we will see 
later, accuracy is achieved only when we use very small dis- 
tances and consequently very large concentrated forces. In 
this way of moment decomposition, the usual behavior of a 
concentrated moment is achieved. For larger distances the 
couple is not an accurate representation of a concentrated 
moment and the results are not so accurate. 
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In the current example we use M = 10 and it is 



decomposed as follows: 

10 = 1000 x 0.01 

a pair of opposite forces of 1000 units acting 0.01 units of 
length apart from each other. In our problem we use 200 
intervals of 0.005 each, one concentrated force is located 
at x = 0.495, and the other opposite force of equal magnitude 
is placed at x = 0.505. So the loading function to be inte- 
grated is defined by 5 partial functions as follows: 



f 1 (x) = 0 


0 < x < 0.495 


J f 2 (x) dx = -1000 


x = 0.495 


O 

II 

X* 

CO 

M-l 


0.495 < x < 0.505 


Jf 4 (x)dx = 1000 


x = 0.505 


f 5 ( x ) = 0 


0.505 < x £ 1.0 



c. Necessary Initial Conditions 

We need to consider now the B.C.'s. For the present 
configuration we have: 
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Shear : 


Y"' = Unknown 

o 


Y " » 

n 


= Unknown 


Moment : 


Y" = Unknown 


Y" 


= 0 




o 


n 




Slope : 


Y' =0 

o 


Y 1 

n 


= Unknown 


Deflection : 


Y =0 


Y 


= 0 



o n 



To start the algorithm we need to determine the initial 
shearing force Y^' and the initial bending moment Y'\ From 
equations (3.13) to (3.16), after replacing the known values 
we get: 

0 = hS 0 + nhY"' + Y" 

2 o o 

Y' = h 2 S. + N 1 h 2 Y M ' + nhY" 
n 1 1 o o 

0 = h 3 S + N h 3 Y + N,h 2 Y" 

o o o 1 o 

This is a set of three simultaneous algebraic equations with 
three unknowns and we solve for the initial conditions. The 
results are: 



Y" 

o 



h(N S„ 

O Z 

nN^ 



nS o ) 



N 



o 



(5.4) 



Y 11 • 

o 




N 1 S 2 



- N 



o 



(5.5) 
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All the terms involved in these expressions are 
defined by the relations (3.13) to (3.22), and they can be 
computed in advance since they are independent of boundary 
conditions. These initial values are only approximations as 
we know, but serve to match the given B.C.'s at the other 
extreme. The results are: 

Shear: Y£' = 11.21219 Exact = 11.25 

Moment: Y” = -1.21219 Exact = -1.25 

o 

Running the program with these initial values, we 
get the results shown on the next page. The deflection appears 
accurate to about 2 percent for the more meaningful values of 
deflection. However, if we reduce the distance between the 
opposite forces, we improve the accuracy. But if we increase 
it, the discrepancies grow enormously even for a very small 
increment, and it does not reflect a concentrated moment 
behavior . 

The slope shows the expected behavior. However, we 
should not expect high accuracy since we have started the 
algorithm with approximate starting points of shear and moment. 
At x = 1.0, we have Y = -0.631. The exact is -0.625. Pre- 
ceding values of slope are expected to be better since here 
we start with an exact known initial point. Similarly for 
the moment we start with an approximation but get more accurate 
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EXAMPLE 2 ‘ 

C LAf^PcC-PQLLE 7 E5A.4 
BOUNDARY CONDITIONS 



Y (0 ) =0.0 Y(1)=0.C 

Y (0)=0.0 M { 1 ) = 0 . C 



X : 


SHEaR 


fCMENT 


SLOPE 


OeFLcCTI JN 


EXACT DEF 


0.00 


11.212190 


- 1.212193 


C.COOOOO 


C. 000000 


0.000000 


0.05 


1 1. 212190 


-C .651584 


-0 .0465 94 


— C.OO 1230 


-0.001327 


0. 10 


11.212190 


-C.C9C975 


—0 .0651 58 


-C. 004190 


-0.00-^374 


0.15 


11. 212190 


C .469634 


-0.C55692 


-C. 007327 


-0.007 734 


0.2C 


11. 212190 


1.020243 


— 0.C13195 


-C. 009290 


-0.009999 


0.25 


1 1.2 12190 


1.590358 


0.C47333 


-C. 008677 


—0 .0U9766 


0. 3 C 


11 .212190 


2.151463 


0 .140890 


-C. 004087 


-0.005624 


0.35 


11.212190 


2.712071 


0.262479 


C. 005332 


0.008828 


0.4C 


11.212190 


3.2 72680 


0.412098 


C. 022630 


C. 020000 


0.45 


11.212190 


2.833291 


0 .589747 


C. 04 756 1 


0.04^297 


0.50 


-988 .73 7500 


-2.106097 


0.7641 77 


C. 081965 


0.078 125 


0.55 


11.212190 


— 5 .045486 


0.504137 


C. 113693 


0.110890 


0 . 6 C 


11.212190 


-4. 48^873 


0 .265877 


C. 132833 


C. 130000 


0.65 


11.212190 


— 2 . 92426 9 


0 .055649 


C. 140755 


0.138359 


0.7C 


11.212190 


-2.363659 


-0.126850 


C. 138867 


0 . 13o 87 5 


0.75 


11.212190 


-2.8C3050 


-0 .230717 


C. i.28570 


0. 126953 


0.3C 


il.212 190 


— 2.2424^0 


-0.406855 


C. 111265 


0.110000 


0.85 


11.212190 


-1 .661830 


-0.504961 


C. 033354 


0.037422 


0.9C 


11 .212190 


-1.121222 


-0.575087 


C. 061239 


0.060625 


0.95 


11.212190 


— C .560ol 3 


-0 .617034 


C. 031319 


0.031016 


1.00 


11 .2 12190 


-0 .000003 


-0 .631098 


C. 000000 


C.OOOCOO 
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as we proceed, since we end up with the expected result of 
0.0 at x = 1.0. 

In the case of the shear force, we see it is con- 
stant throughout the beam length as expected but it is shifted 
by an amount of 0.0378 from the exact. It is interesting to 
note that at the point of the loading we have a very strong 
shearing force. For a pure moment acting at this place, we 
shall neglect this result. However, if we do have an actual 
couple than a big shearing at this point should be expected. 
Note that the algebraic sum of shear forces gives: 

11.21219 - (-988.7875) = 1000.0 



which is indeed the couple acting at this point. 

5 . Example 3: Pin-Roller Beam 

Refer to Fig. 5.1c. This is a statically determinate 
simple supported beam subjected to a triangular type of load 
distributed over the central portion of the beam as indicated 
in the figure. We need to determine the shear, moment, slope 
and deflection of the beam. 

6 . Solution 

First of all, some comments about this configuration 
are in order. This example was chosen in order to emphasize 
the ease with which we can switch from one kind of load dis- 
tribution to another. As explained at the beginning of this 
work, it is the fact that the algorithm solves a "new" differ- 
ential equation at every single step that allows us to deal 
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with different loads. Therefore it is possible to have a 
different kind of loading per subdivision. Another reason for 
selecting this problem is that here we deal with linear forcing 
functions and we are going to perform four successive integra- 
tions of a linear function. This will lead to a polynomial 
of fifth degree and a bigger cumulative error is to be expected. 

Proceding with the solution, the distributed load 
f(x) to be integrated is given by four different functions: 

f , (x) = 0 0 £ x < 0 . 25 

f 2 <x) = 4x-l 0.25 £ x < 0.50 

f^(x) = -4x+3 0.50 £ x < 0.75 

f^(x) = 0 0.75 £ x £ 1.0 

As indicated above, the computer program is such that 
it can switch from one kind of loading to the other at the 
specified point. Now, in order to start the algorithm we 



need to 


determine the 


initial 


points from 


the B 


.c 


. ' s of a 


simple 


supported beam. 


These 


are : 










Shear: 


Y "» — 


Unknown 


Y 1,1 




Unknown 






o 




n 








Moment : 


Y" = 


0 


Y" 


= 


0 






o 




n 








Slope : 


Y 1 = 

o 


Unknown 


Y 1 

n 


= 


Unknown 




Deflection : 


Y 

o 


0 


Y 

n 


= 


0 
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Now, from the known relationships after we substitute in the 
given values, we get: 



0 = hS~ + nhY'" 

2 o 

0 = h 3 S + N h 3 Y"' + nhY ' 

o o o o 



This is a pair of simultaneous algebraic equations in two 

unknowns. We solve for Y"' and Y'. The results are: 

o o 



Y ' = (£) 2 (N S 9 - nS ) 

o no Z o 



(5.6) 



Y ** 1 
o 




n 



(5.7) 



Similarly, the computer supplies the values of the 
variables involved in these relations and these are as follows 



Shear: Y"' = -0.12499 Exact = -0.125 

o 

Moment: Y' = 0.014972 Exact = 0.014974 

o 

which shows a remarkable accuracy. Running the program with 
these initial values we obtain the results shown on the next 
page. A comparison between the approximate deflection and the 
exact reveals an error of about 2 percent at the point where 
the largest deflection takes place. In this example, we 
expected a larger error, but the algorithm appears to perform 
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EXAMPLE 3 ‘ 

P IN— ROLLER 3E*M 



BOUNDARY 


CONDITIONS : 


Y (0) =0. 


0 Y 1 1 ) = 0 


.0 








M ( Q )=0. 


0 M ( 1 ) = 0 


.0 




X : 


SHEAR 


NO, ME NT 


SLOPE 


DEFLECTION 


EXACT DEF 


0.00 


-0. 124599 


C.COOOOO 


0 .C14972 


C. 000000 


0.000000 


0.05 


-0.124599 


-C.0Co250 


0 .C14316 


C. 000746 


0.900746 


0.10 


-0.124999 


-0.012500 


0 . C 14347 


C. 00 1476 


0.001477 


0.15 


-0. 12 459 9 


-C. 016750 


0.C13566 


C. 002175 


0.002176 


0 . 2c 


-0.124999 


-0.025000 


0 .C 12472 


C. 002 32 7 


0.002 3 23 


0.25 


-0. 124599 


-0.021250 


0.C11066 


C. 00341 7 


0.003418 


0.30 


-0. 115599 


-C. 037415 


0 .C093*3 


C. 003929 


0.003920 


0.35 


-0.104999 


-0 .04^030 


0 .007333 


C. 004347 


0 . 00*322 


0.40 


-0. C75999 


-C.C4 7745 


0.CQ5O57 


C. 00465 d 


0.004609 


0.45 


— 0 . C -*4999 


-C. 050910 


0 .C02534 


C. 004349 


0 .004 76 6 


0.50 


0. COCCOO 


-0.052075 


O.COOOOO 


C. 004914 


0.00*735 


0.55 


0. C45CU0 


-C. 050910 


-0.C02534 


C. 004349 


0.0047o6 


O.o C 


0. 33CCOO 


-C. 047745 


-0 .C05057 


C. 004653 


0.0046C9 


O.c.5 


0.105000 


-C. 043030 


-0.C07333 


C. 004347 


0.004 322 


0.7C 


0. 12CCOO 


-C. 037415 


-0 .009346 


C. 003929 


0.003920 


0.75 


0.125000 


-0.031250 


-0 .Cl 1066 


C. 003417 


C. 0034 18 


0 . 3 C 


0. 125000 


-0. 02 5 00 0 


-0 . C 124 72 


C. 002327 


0 . 00 *1 3 l . 8 


0.85 


0. 125GOO 


-0.013750 


-0 .013505 


C. 00217;. 


0.002176 


0.90 


0 . 125000 


-0.012500 


-0 .0 14347 


C. 001476 


0.001477 


0.95 


0. 125C00 


-C.CCO250 


-0 .014316 


C . 00 0 7*6 


0 .000 7 4o 


1.00 


0.125000 


C.COOOOO 


-0 .014972 


- C. 000000 


0.000000 
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fairly well even in this case. Nevertheless, the error should 
be expected at the "free" (unknown) B.C.'s since the other 
extremes are fixed. Here we make use of the fact that the 
loading is symmetrical and the results should be symmetrical 
too. This is actually the case in the computer printout, so 
we can conclude that the results are indeed accurate. 

7 . Example 4 : Clamped-Clamped Beam 

Refer to Fig. 5. Id. This is a statically indeterminate 
beam clamped at both ends and loaded with a totally arbitrary 
discontinuous load known only as a table of values at discrete 
points as shown on the next page. We shall determine the shear, 
moment, slope and deflection of the beam. 

8 . Solution 

Before we proceed with this example, let us remind 
ourselves that we are not restricted to handle only exact- 
integrable forcing functions. We can deal with arbitrary 
functions as well, and this example intends to illustrate the 
procedure. In this problem, the computer program reads in 
data from a table of values instead of obtaining them by 
evaluating a given function. With this data the program per- 
forms trapezoidal integration in the usual way. The forcing 
function is defined by three other functions as: 

f x ( x ) =0 0 < x < 0.3 

f 2 (x) = (from data deck) 0 . 3 _< x £ 0 . 7 

f 3 (x) = 0 0.7 < x < 1.0 
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OAT 5 C2CK FCF 


2X.i'2PLi 4 






OISTANC'J 


LOAD 






o. :c 


0 . 300333 


U . 5 6 


6 . 4000o0 


o. : l 


C. 000030 


0.57 


'6 . 2'JoOOO 


0.32 


0 .300030 


0.5c 


7 • b J 0 0 u 0 


0. J3 


C .030000 


0.59 


6 .9 JOOOO 


o. :4 


3 • 0030 l/0 


3.30 


6 . 300000 


3.0 0 


0 . 000030 


0.51 


4. 1000 00 


J » v 6 


0 .0503 00 


0.62 


4.000000 


0.37 


C. 030030 


0.52 


4 .600000 


0. 38 


0 .000300 


0.64 


3. 503000 


0.39 


0 . JOOoOO 


0.65 


4.000000 


J . 1 C 


G . 030030 


0.56 


3 .200000 


0 . L 1 


0 • 0000 00 


0 • o 7 


3.500000 


0.12 


C .030000 


0 • o 3 


3.200000 


3.10 


C. 3000G0 


0.59 


3 .4000 JO 


3.14 


0 .000330 


0. 7C 


3. JOOOOO 


3.10 


0 . J00GO0 


0 . / i 


C. 0000 JO 


0. 16 


3 . Go 0030 


0 . / 2 


0 .000000 


0.17 


0. 000330 


0. / 3 


C.OJOOOO 


o . i e * 


0 .0003 o J 


0.74 


c. 000000 


0. 1 9 


0 . u'OOOOO 


0.7 5 


0 . JOuOGO 


3. 2C 


G .300000 


0. It 


G. 000000 


0.01 


0 • 3 3 0 0 5 0 


0. / 7 


C.OJOOOO 


0.22 


G. 00 0000 


0 . 7£ 


c .000 0 00 


0. 2 J 


0.0030 JO 


0.7 9 


0. oOOOoO 


0.2 4 


0 .u 0 J jOO 


0.30 


C.cJ^OOJ 


3.2 0 


G . J JOGoO 


u. 2 1 


c .000000 


3.2c 


0 . )303 JO 


0.3 2 


C. JOOOOO 


3.27 


0 .o^oOoO 


0.23 


0.000000 


0.26 


G. . o 0 o 3 0 


0.34 


0.0 JOOOO 


0. 29 


0 . 303030 


0.35 


C . 0 300 00 


3. 2C 


l.OOOOoO 


0.3c 


0 .OOoOOO 


0.2 1 


l.oOOOOO 


0.37 


C.OJOOOO 


0 . 3<l 


1 .3300 30 


0.38 


C. JOOoOO 


3.20 


2 . oOuO JO 


0.39 


G .030000 


0.24 


2.200033 


0. J C 


C. JO JO 00 


0. 35 


2 . o 0 5 0 0 0 


0 . J 1 


C.OJOOOO 


3. J 6 


3.2 JooOO 


0.32 


0 . OOOOOO 


0.37 


3 .000003 


0.5 3 


C. J JOOOO 


0.3c 


3 .400030 


0. 34 


0.000000 


3.39 


3 . 6 3 u o 3 0 


J . J 5 


3 . JJUJUO 


3 . 1 C 


o . 030 0 30 


0 . j 6 


G • J JOOoO 


0.71 


5 .3 JOOOO 


J • 3 7 


0 . OOoOOO 


0.4^ 


6.2000^0 


0.38 


0 . OOOUuO 


0.75 


6 . 600000 


0.3 9 


C. OOOOOO 


3.74 


6 .700300 


1 . JC 


0 .uOoGOO 


0.7 5 


6 . 4003 00 






0 . 4 fc 


6. 200000 






3.77 


7 .3000 oO 






0.7c 


6 . 900000 






3 . ♦ 9 


6.300000 






0. 5C 


6. 1030 JO 






0.5 1 


6.030000 






0.52 


5.330000 






0.5 5 


6.200330 






0.54 


6 . 6 OOo 00 






0.55 


6.700000 
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Unfortunately, no exact solution has been obtained for 
this case and we are not able to compare the computer results 
with the exact. Nevertheless some conclusions can be drawn. 
The B.C.'s for a clamped-clamped beam are as follows: 



Shear: Y"' 

o 



Moment: Y" 

o 



Slope: Y^ 



Deflection: Y 

o 



Unknown 

Unknown 

0 

0 



Y M » 

n 



Y " 



n 



Y 



n 



Y 



n 



Unknown 

Unknown 

0 

0 



To determine the initial conditions, we use the known 
relationships and substituting in the known values, we obtain: 

0 = h 2 S. + N-, h 2 Y + nhY" 

1 1 o o 

0 = h 3 S + N h 3 Y'" + N,h 2 Y" 

o o o 1 o 

This is a pair of simultaneous algebraic equations 
with two unknowns. Solving for the initial conditions we 
get: 



Y" 

o 



n (N S 0 - nS ) 

Q Z Q 

nN, - N 
1 o 



(5.8) 



Y n i 

o 



S - 

o 

nN 1 



N 1 S 2 



- N 



(5.9) 
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For this example the values supplied by the computer 



are : 



Shear: Y’" = -0.948411 

o 

Moment: Y" = 0.230507 

o 

The computer printout is shown on the next page. Here 
we have another column to show the arbitrary loading. Note 
that we are using 100 intervals, but actually every fifth is 
printed . 

Perhaps the only way to analyze these results would be 
to check if the B.C.'s have been met or not. The results show 
that they have. Another clue to assure some accuracy would be 
to look at the maximum deflection and the minimum slope. The 
loading was intentionally concentrated on the middle of the 
beam so we can expect the highest deflection and minimum slope 
in this neighborhood since we have similar B.C.'s at both 
ends. We can say that this too checks. Consequently we see 
that the results are likely to be trusted. 
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EXAMPLE 4 

CLA^PEC-CLAMP EC 3 EAN 



BCOAC AF Y 


CGM C I T I * ,'n: S : 


Y ( u ) - J .0 


Y ( 11=0.0 










Y 1 ( 0 ) = 0 . J 


Y' ( 11=0.0 






X : 


l:pc 


it- TAR 


f'G.’l - IT 


slop: 


DEFLECTION 


0. JC 


0 . ccccoo 


-C.94C41 1 


0 .230507 


C. OGOOUO 


0 .000000 


0. 05 


0.0 JCCJO 


-C. 948411 


0 . 1S3J36 


C . 01 0340 


0.000268 


0. 1C 


0. CCCCJU 


-C.94c4i* 


0 . 1 J 5o66 


C. Cl 3209 


0.000994 


0 • 1 5 


0. CCCCOO 


-C. 942411 


0 .033245 


C. 023906 


0.002059 


0 • 2 C 


o.occoo 


-0 . 94 c4 1 1 


0 .040325 


C • 0<i 71 j3 


0.005344 


0.25 


C. CCCCOO. 


— 0.948 4 i 1 


-J.CG6596 


C. 027939 


0.004732 


U.3C 


1. 6'wuCOO 


-C.94j411 


-0 .05 3901 


C. 026474 


0.006 103 


0. -> 5 


Z.OCCOuO 


. 84 1 


-0 . C940o7 


C. 02262 ‘J 


0. J07339 


0. 4C 


o • w^COjO 


-u. 67 2^1 i 


-0.13 74 a 7 


C • oi 6630 


0 .008330 


0.45 


6. 4 - C G 0 0 


-0.35 / 4 L i 


-0 . lbi »28 


C. 009094 


0. J 0 o 9 7 9 


0 • dQ 


0 . i J oooo 


— c • o 3 u 9 1 1 


~ 0 .173116 


C. j006l4 


0.00 *224 


0. 35 


6.7ICC 00 


C. 27 9 086 


-0 .16 7032 


-C. 007952 


0.009039 


0. oC 


6 «JuCjuO 


C • 605 08 5 


-0 . 14 4 0 53 


-C. 01531 4 


0. J08H40 


0.65 


4 . :c ccoo 


o • c 2 5 o J j 


— 0 . 1033 41 


-C. 022^07 


0.007^83 


0.7C 


3. u w L L O J 


C .996 580 


-0 .06 3138 


-C. 026542 


0.006255 


0. 75 


u. w j C C 0 0 


i. 01 1530 


-0 .012634 


-C. 02344 0 


0.004370 


0. JC 


o. c : cc ou 


±•01*530 


0 .037344 


-C. 027810 


0. 005454 


0 . o 5 


J»gjI»L>UG 


1 .011580 


0 .033472 


-0.024651 


0.002 132 


C . 9 c 


0. OJOoOO 


i.011530 


0.139050 


- C. 01 39o 3 


0.001032 


0.95 


0. w 2 0 00 0 


1 .011530 


0. 139629 


- C . 01 0 746 


0.000279 


L .00 


0.0 JOCOO 


1 • 0 L 15oO 


0 .240207 


C. 000000 


-C. 000000 
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C. THE COMPUTER PROGRAM 



file: SLOPE WATFIV A 1 



SJCE XKEF 

THE COMPUTER FRCC.RAM 

C«* $***❖❖$$❖$ 

C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

MAIN PROGRAM ❖$$=£$ 



this program calculates the shearing fcrce. bencing MOMENT, slope and 
DEFLECTION of eeams by performing four succesive integrations gn 
THE LOADING FLNCTIQN. THE FOLLOWING CONFIGURATIONS ARE CONSIDERED: 

1 . -CANTILEVER (CLA MPEC-FRES ) BEAM 

2. -CLAMPE0-PQLLER EE A M 

3 • — S I MPLc SUPCRTEO (PIN-ROLLER) EE A M 
A, —CLAMPED— CLAMPED BEAM 

The LOADING CAN BE ANY number OF DISTRUBUTED ANC/CR CONCENTRATED 
FORCES. CONCENTRATED MOMENTS CAN BE DECCMPCSEC INTO A CCUPLE AT THE 
MINIMUM POSSIELE CISTANCE ACTING AT THE PCINT OF APPLICATION CF THE 
MOMENT. 

THE PROGRAM CONSISTS OF THREE ROUTINES: 

1. SUBROUTINE SLOPE 

2. SUBROUTINE INCCM 

3. SUBROUTINE LCAC 

HH ICH ARE DEFINED IN THE FOLLOWING PARAGRAPHS 



C READS IN NUMBER CF STEPS CE FINED BY THE INTEGER VARIALBLE N 
C CALLS SUBROUTINE INCCM FOR INITIAL CONDITIONS 
C CALLS SUBROUTINE SLCPfc FOR FINAL RESULTS 
C 

SUBROLTINE INCCM $** 2 **^*^ ^ 

c 

READS IN THE INTEGER VARIABLE K WHICH SPECIFIES THE KIND OF PROBLEM 
WITH THE FOLLOWING CODE: 

K=I CANTILEVER BEAN (INITIAL CONDITIONS MUST BE SPECIFIEC) 

K = 2 CLAMPED ROLLER BEAM 
<=3 C LA MPED— CLAMPED BEAM 
K=4 SIMPLE SLPPCRTEC BEAM 



C 

c 
c 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

SUHRCUTINE SLOFE 

C PERFORM TwO TASKS: l. CALCULATES ANO SUPPLY INCCM WITH THE NECESARY 

PARAMETERS TC DETERMINE INITIAL CCNCITICNS 
2. EXECUTES THE SOLUTION CF THE PROBLEM 



YO 

Y 10 

Y <c C 

Y 30 



THE PARAMENTERS ARE DEFINED AS FOLLOWS: 
CALL tNC0N(N.YC.YlC.Y20 .Y3C) 

N = NUMBER CF STEPS 

INITAL CEFLECTICN 
INITIAL SLOPE 
INITIAL BENDING MOMENT 
INITIAL SHERING FCRCE 



THIS SUBROUTINE SUPPLIES THE INITIAL CONDITIONS BY CALLING SLGPE. 

fcr a cantilever eeam this statement is nct executed 



this subroutine is a the coding of the relationships given by 

EOLATIONS 3.13 TO 3.32. THE SLMMATICN TERMS ARE CONTAINED IN 

the common statement 

THE CALLING PARAMETERS ARE DEFINED AS FOLLOWS: 

CALL SLOPE ( N. YO ,Y IC .Y23 . Y3C. J ) 

N = SAME AS IN INCCM 

YC = " 

Y1C = " 

Y2C = ** 

Y3C = » 

J = TAKE two values: j=o: tc find initial conditions only 

J= I : TC EXECUTE THE SOLUTION 



SUBROUTINE LOAD 

c 

C SUPPLIES SLOPE WITH THE NECESARY INFORMATION OBTAINED FRCN THE 
C LOADING FUNCTION. IT CAN EE AN ACTUAL FUNCTION CR A DATA DECK. 
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file: slope 



fcATF |V 



A | 



MUST dE WRITTEN SPECIFICALLY FOR EACH NEW PRCfi LE M . «HEN OEAlING 
with CONCENTPATcO LOACS. the VARIABLE PINT in RCUTInE SLOPE is 
IS SET EQUAL to THE CONCENTRATED LCAO DIRECTLY • AT THE POINT OF 
ACTION. 

The CALLING parameters ARE CEFINED as FOLLOWS 
Call LOAO( e • H 2 ) 

CURRENT VALLE OF THE VARIABLE X (LENGHT) GIVEN BY ROUTINE 
H2= CONTAINS THE VALUE OF THE LOADING FUNCTION AT X = 0 



C MAIN PROGRAM 



REAL YO.Y lC.Y2C.YJQ 
INTEGER J.N 

COMMON SO. SI . 32. SNO.SN1 .RN.H 
REAO( 5,10)N 
10 FORMAT! | 3 ) 

H=| ./FLJA T(N- 1 ) 

RN=FL JA T ( N- l ) 

CALL INCOMN . YO. YiO, Y20 . Y30 ) 
J=0 



CALL SLOPE(N. YO. Y10.Y20. Y^O.J ) 

STOP 

END 



C 

SUBROUTINE LCA0(e.H2) 



slope 



REAL a.H2 
FIX)— .... 

H2=F( E ) 

RETURN 

ENO 

C 

cs S* ^ 4: ££ $ * $ £ £ S 

SUB ROUT I NE INCCN (N.YO.YIC.Y 20 • Y 3 C ) 

REAL Y0.YlC.Y20.Y20 
INTEGER J.K.N 

COMMON S0.Sl.S2.SN0.SM .RN.H 
REAU(5.lO)K,YC. Y10.Y2C.Y3C 
10 FORMAT! I | .4F 10. A ) 

IF! !K.EQ. C ).QR. ( K.GT.4) ) GO TO 60 
IF! K. EQ. I ) GO TC SO 
J = l 

CALL SLOPE ! N . YC . Y l 0. Y20 . Y30 . J ) 

IF( K-3 ) 2 0 .30.40 

c clamped-rcller 

c ^ 

20 WRI TE( 6. 2E ) 

2S FORMA T( • l •.• EXAMPLE 2 •// l X. *C LAMPcD-RCLLER EEA M* // I X . ^BOUNDARY CON 
COITIONS: Y(0)=0.0 Y(| ) = 0.0 *//27x , *Y (0>=0.0 M ! I ) =0 . o • ✓/ ) 

Y2Q=H^( S2^SN0-RN^S0 )/ (RN^SNl-SNO) 

Y30- ( S0-SNl^S2)/( RN«SM-SNO ) 

RETURN 

C clamped-clamped 

30 WRI TE ( 6 * 3 S ) 

2S FORMA T !• I •.• EXAMPLE 4 •// I X. *C LAMPEC-CLAMPED 0E AM • / / l X , • B OL NO AR Y CO 
SNOITIONS: Y!0)=0.0 y! I ) =0 . 0 • //2 sx , • Y !0) = 0.0 Y ! I ) =0 .0 •✓✓ ) 

Y20=H4!SO^SNI-Sl^SNO) / ( RN^SNO -SNt^SNt ) 

Y JO — ( SI ^SN l -RNiSO )/(RN*$NO-SNl*SNl ) 

RETURN 

C PIN-ROLLER 

C S ^ ^ S *=£:**: £3: £ S ^ S - 3: S « S S ^ - 
40 WRI TE( 6. 4 E ) 

4 s form a T( • i • • 'Example 3 • // 1 x. •«* i n-roller beam* // ix . *bcunoary conoiti 
5 0NS: Y ( 0 ) =U . 0 Y ( l ) = (J.O'//27x. *M( 0)=0.0 M(l)=0.0*//) 

YiO=( (H/RN 2 ) - ( SN0-S2-RN-S0 ) 
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file: slope 



• ATF I V 



A 



Y30 = -S^/fiN 
«E T UflN 

C CLAMPED-FPEE 

5 C «RI TE ( 6, 55 ) YO* Y 1 C . Y2C • Y2C 

55 FORMAT^ 1 ‘.'EXAMPLE I •// IX . »C LAHPEO-FREE BEAM* // IX « * BCUNOAP Y CON 
COITIONS: Y(0)='.F8.4,4x, »Y (0)=* ,F8,4,4X,//27X, ’MIOls'.Fa.A, 

*4X, •V(O)s».F0.4//) 

GO TO 90 

60 MR I TF ( 6 • 65 ) 

65 F 09 MAT( 1 X » *N FROM I TC 4 ONLY • / ) 

90 RETURN 
ENO 
C 

c ^ ^ £ $■ =* ■$ S =}=£ ^ 

subroutine slope (N.Y0.Y10.Y 20 ♦ YJC« J ) 

REAL SUO 0 . SUO 1 • SU02.SU03. SU 1 1 . SU 1 2 • SU I 3 . SU 22 . S U2 3 . SU33 • R I N T , 

*RINT1 •RINT2.RINT3,Hl,h2.RIta.SNMl*SNM2 
INTEGER I • J • K .L.N t M,NSlK«N SOL .NSlU.IMl 
COMMON S0.S1,S2.SN0,SM «RN«H 
SUOO=G. 

SUO 1=0 . 

S UO 2 = 0 . 

SOO 3=0. 

SU1 1 =0. 

SUI 2-0. 

SU1 3=0. 

SU2 2 = 0 . 

SU2 3 = 0. 

SU33=0. 

R IN T = 0 • 

RINT1 =0 . 

R IN T2 =0 • 

RINT3=0. 

NS 1 K=0 
NSO L = 0 
NS1 L=0 
H I = 0.0 
H2 = 0 • 0 

IF(J.GT.O) GO TO 100 
WR1 TE(6,2CJ 
100 00 400 i = 1 • N 

1 Ml = 1 - 1 
RI=FLOAT( IM1 J 
9=R I*H 

CALL LOAD ( 0 »H2 ) 

1F1 I .EG. 1 )GC TO 200 
RINT=0.5~FS<H1>H2) 

200 RINT4=RINT3 
RINT3=R1 NT2 
R IN T2 = RI NT 1 
RINT1 =PI NT 
C 

SU33=SU334P INTI 
S3 = SU 33 
C 

SU23 = SU23 4PI NT2 
SU2 2 = SU2 2 +SU23 
S2= . 5-SU33> SU22 
C 

K = I -2 

I F ( I.LE. 2 ) K = 0 

SUI 3= SU 1 3 +RINT3 

SUI 2= SU 1 2 +SL 1 3 

SUI 1 =SU1 1 *S01 2 

NSl K = N S 1 K ♦ K 

SNM 1=FL0 AT (NS IK J 

SN1 =0.5^FI*SNM1 

Sl= . 25FSU22+SL22+SU 1 1 
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FILE 



SLOPE 



HATF IV 



A l 



L=I-3 

I F ( I.LE. 3 ) L“0 
SUO J=SU0 3 4R IM4 
SU0 2=SUQ 24SU03 
SUO 1 = SUG 1 4 SUQ 2 
SUOO^SUGO 4SUO 1 
NSIL=NS1L4L 
NSOL=NSOL 4NSI L 
SNM2 = FL3A KNSOL ) 

SNO =0.25iRl+SNMlfSNH2 
SO= • I 2S4:SU334 .7SSSU22 + 1 .5$SU1 l+SUQO 
I F ( J • GT • 0 ) GO TQ 300 
C 

Y3=S3>Y30 

Y2-H* ( S2 + R I«Y30) 4Y2C 

Yl - ( ( SI 4SN 1- Y3C) + I £ Y 2 C )-H4 Yt 0 
Y=( ( (S0 + SN0-Y3Q) *H+SN 1 * Y 2 C ) -H ♦ » I * Y 1 O ) 4 YC 
*RITE{6, lC)b.Y3,Y2. Yl . Y 
300 hi=h2 
400 CONTINUE 

i f ( j • gt • o ) go to ggq 

WRI TE (6.95) 

10 FORMAT! 1X.F4.2.2X.4F12.6/) 

2 0 FOR MAT( 2X • • XI * . 9X . • SHE AR • , bX . 'MOMENT ' .6X. • SLOPE* . 5 X • • OEF L EC T I ON • 
RR FORMAT! *| •) 
qg<5 return 
end 
c 

SENTRY 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



Based upon the research carried out in this thesis and 
the results obtained, the following conclusions can be drawn: 

1) First of all, the method, in the way it has been 
developed, shows the fact that the solutions are totally 
independent of each other. Any value of the unknown can be 
found without generating all previous ones. A closer look at 
equations (3.13) to (3.16) reveals that we can apply any of 
these relationships at any point i (0 < i < n) directly, pro- 
vided we know the initial conditions and the summation terms 
which can be generated in advance. In all of these equations 
there is only one summation term which depends on the given 
function f(x), the other summations are series of integers 
totally independent of the given problem. This fact represents 
a good saving in terms of computer time if it is conveniently 
exploited . 

2) As it has been shown in the examples, the power of 
the method is perhaps its ability to deal with arbitrary func- 
tions and this is important since many engineering problems 
lead to these kinds of functions. 

3) So far, there has not been any error correction in a 
strict sense. The only corrective measure has been step size 
reduction. However, equation (4.9) shows that correction can 
be performed provided we know how the local error behaves. 
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When integrating constant and/or linear functions, (4.9) 
applies, but for other cases it does not. In any case, equation 
(4.3) indicates that the local error depends only on the second 
integral which in turn is one degree higher than the given 
function itself. This suggests the idea that if we know how 
the original function behaves, linearly, quadratically , etc., 
we could, to some extent, predict the error behavior and carry 
out a correction. Further research is recommended in this 
particular case. 

Nevertheless, as shown in the examples, some accuracy has 
been achieved even working in single precision. For more 
complicated problems involving complex functions and requir- 
ing very accurate solutions, double precision is still a good 
possibility . 
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